function [z1,z2,z3,z4] = plant_SPDE2D(M,L,c,c2,sigma,time,dt,z10,z20,z30,z40)
if nargin == 7
    z10 = zeros(M+1,M+1);
    z20 = z10;
    z30 = z10;
    z40 = z10;
end

assert(size(z10,1)==M+1 && size(z10,2)==M+1);
assert(size(z20,1)==M+1 && size(z20,1)==M+1);
assert(size(z30,1)==M+1 && size(z30,1)==M+1);
assert(size(z40,1)==M+1 && size(z40,1)==M+1);

z1 = cell(length(time),1);
z2 = cell(length(time),1);
z3 = cell(length(time),1);
z4 = cell(length(time),1);
%%
t = 0;
lambda = zeros(M+1,M+1);
for m = 0:M
    for n = 0:M
        lambda(m+1,n+1) = -4*c2*(m^2+n^2);
    end
end
sqrt_dt = sqrt(dt);

if(time(1)==0)
    z1{1} = z10;
    z2{1} = z20;
    z3{1} = z30;
    z4{1} = z40;
    i = 2;
else
    i = 1;
end

z1_state = zeros(M+1,M+1);
z2_state = z1_state;
z3_state = z1_state;
z4_state = z1_state;

c_200 = c/L^2;
while t <= time(end)
    for m = 0:M
        for n = 0:M
            if m > 0 && n > 0
                z1_state(m+1,n+1) = z1_state(m+1,n+1)+lambda(m+1,n+1)*z1_state(m+1,n+1)*dt+randn()*sqrt_dt*sigma;
            end
            
            if m == 0 && n==0
                z2_state(m+1,n+1) = z2_state(m+1,n+1)+c_200*dt+randn()*sqrt_dt*sigma;
            else
                z2_state(m+1,n+1) = z2_state(m+1,n+1)+lambda(m+1,n+1)*z2_state(m+1,n+1)*dt+randn()*sqrt_dt*sigma;
            end
            
            if m > 0
                z3_state(m+1,n+1) = z3_state(m+1,n+1)+lambda(m+1,n+1)*z3_state(m+1,n+1)*dt+randn()*sqrt_dt*sigma;
            end
            
            if n > 0
                z4_state(m+1,n+1) = z4_state(m+1,n+1)+lambda(m+1,n+1)*z4_state(m+1,n+1)*dt+randn()*sqrt_dt*sigma;
            end
        end
        
    end
    t = t+dt;
    if(t>=time(i))
        z1{i} = z1_state;
        z2{i} = z2_state;
        z3{i} = z3_state;
        z4{i} = z4_state;
        i = i+1;
    end
end